function corrected_freq = max_spectrum_correction(signal, fs, de, max_iter)
% 最大值估计法频谱校正函数
% 参数:
%   signal   - 输入的中频信号 (时域)
%   fs       - 采样频率 (Hz)
%   de       - 允许的最大频率误差阈值 (Hz)
%   max_iter - 最大迭代次数 (默认100)
% 返回:
%   corrected_freq - 校正后的频率估计值 (Hz)

if nargin < 4
    max_iter = 100; % 默认最大迭代次数
end

n = length(signal);
fft_result = fft(signal);
freq_res = fs/n;
freq_axis = (0:n-1)*fs/n; % 生成正频率轴

% 转换为双边谱
if mod(n,2) == 0
    freq_axis(n/2+1:end) = freq_axis(n/2+1:end) - fs;
else
    freq_axis((n+1)/2+1:end) = freq_axis((n+1)/2+1:end) - fs;
end

magnitude = abs(fft_result);

% 找到主瓣内的最大值和次大值谱线
[~, max_idx] = max(magnitude);

% 寻找相邻次大值 (确保在主瓣内)
if max_idx == 1
    second_max_idx = 2;
elseif max_idx == n
    second_max_idx = n-1;
else
    left_neighbor = magnitude(max_idx-1);
    right_neighbor = magnitude(max_idx+1);
    if left_neighbor > right_neighbor
        second_max_idx = max_idx - 1;
    else
        second_max_idx = max_idx + 1;
    end
end

% 初始化参数
ft = freq_axis(max_idx);
fx = freq_axis(second_max_idx);
k1 = magnitude(max_idx);
k2 = magnitude(second_max_idx);

% 迭代修正过程
for iter = 1:max_iter
    m = (k2 - k1) / (k1 + eps); % 防止除以零
    denominator = m - 2;
    
    if abs(denominator) < eps
        break;
    end
    
    e = ((m-1) * (fx - ft)) / denominator;
    
    if abs(e) <= de
        break;
    end
    
    % 调整谱线位置
    delta = sign(e) * de;
    ft = ft + delta;
    fx = fx + delta;
end

corrected_freq = ft;
end